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Abstract The CESAM code is a consistent set of pro- 
grams and routines which perform calculations of ID 
quasi-hydrostatic stellar evolution including microscopic 
diffusion of chemical species and diffusion of angular mo- 
mentum. The solution of the quasi-static equilibrium is 
performed by a collocation method based on piecewise 
polynomials approximations projected on a B-spline ba- 
sis; that allows stable and robust calculations, and the 
exact restitution of the solution, not only at grid points, 
even for the discontinuous variables. Other advantages 
are the monitoring by only one parameter of the accu- 
racy and its improvement by super-convergence. An au- 
tomatic mesh refinement has been designed for adjusting 
the localisations of grid points according to the changes 
of unknowns. For standard models, the evolution of the 
chemical composition is solved by stiffly stable schemes 
of orders up to four; in the convection zones mixing and 
evolution of chemical are simultaneous. The solution of 
the diffusion equation employs the Galerkin finite ele- 
ments scheme; the mixing of chemicals is then performed 
by a strong turbulent diffusion. A precise restoration of 
the atmosphere is allowed for. 
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5 Place 



1 Introduction to CESAM. 

Within the limitations due to electronic degeneracy, CE- 
SAM allows the computation of the quasi-static evolution 
of stellar models as long as the assumption of quasi-static 
equilibrium remains valid, that is to say, until the ex- 
haustion of oxygen in the core. The modular structure 
of CESAM facilitates the choices among several physical 
formalisms, for equation of state (hereafter EOS), convec- 
tion, opacities, diffusion coefficients, etc... Many nuclear 
networks and initial mixtures are available which allow to 
optimise the physical description according to the kind 
of model and evolutionary phase of interest. Mass loss 
and infall of planetoids are also implemented. 



The available packages. Earlier versions, CESAM2-3-4, 
were programmed in F77 and CESAmS in F90. Though 
obsolete, cesam4 and cesamS are available at: 
|http : //www . obs-nice . fr/morel/CESAM| 
In the early 2000's, CESAM was re-programmed in F95 
and named CESAm2k. Three versions are now available: 

— The fixed version]^ and the fixed COROT version" , 
limited to 3a burning, available respectively at: 

,http : //www ■ obs- nice . f r/ cesam/ 
|http : / /per so . obspm . fx/ ~lebreton/Modeles 7CESAM . html | 

— The "/3 versiori^ more complete, not fixed, still in de- 
velopment (not free of bugs) is hereafter signalled by 
the flag {(3). It includes evolution up to oxygen burn- 
ing, diffusion of the angular momentum and other 
developments of minor importance. It is available at: 

|http : //www ■ obs-nice . fr /morel/ CESAM | 

Each package contains five directories and two files: 

— "SOURCE" , contains the FORTRAN sources. 

— "EXPLOIT", contains programs to exploit the models 
and examples of input files. 

— "SUN_STARJDATA" , contains physical data and pro- 
grams for their implementation. 

— "TESTS" , contains programs performing various checks. 
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— "SCRIPTS", contains scripts to be used for the im- 
plantation and operating, and a MAKEFILE. 

— "aide_niem2k . ps" , a short guide of directions for use, 
the "aide-memoire" (hereafter Paper 2). 

— "cesam2k.ps", a complete description of the numer- 
ical aspects and physics implemented, the "notice" 
(hereafter Paper 3). 

The source is structured in 14 modules: numerical rou- 
tines, opacities, convection, etc. All the F95 routines of 
the source are compiled once. The requirements for the 
calculations are read in external files to be supplied by 
the user. The run is interactive. Messages displayed in 
French, or in English, allow the control of calculations. 
The extent of convection zones, a H~R diagram, and the 
profiles of temperature, pressure, luminosity and abun- 
dances are displayecd on line. CESAM has been especially 
designed to facilitate the implementation of various phys- 
ical constants, opacities, EOS, atmosphere, nuclear net- 
works etc. So, its overall structure is separated in two 
spaces: 

1. A "physical space" where the coefficients of the dif- 
ferential equations are written in a form close to their 
physical formalism. 

2. A "nwrnerica/ space" where the differential equations 
are formally solved. 

Therefore CESAM allows to implement physical processes, 
and physical data, without any knowledge of numerical 
methods involved for the solution of the equations. For 
the physics, the use of generic routines makes the reading 
of the algorithms easier. 

Units and values of physical constants. CESAM uses cgs 
units except for the mass, radius and luminosity ex- 
pressed in solar units. Two sets of fundamental constants 
are im plemented which correspond to widely used values 
(IClaytonl.ll968l:IChristensen-Dalsgaar(lll988l : lLi et all . 
119941: ICoxl . l2'000ir For each calculation, one of these sets 
is chosen as the unique source of fundamental constants. 
Other constants are initialised locally, e.g. the mass ex- 
cesses, in the routine performing the calculations of ther- 
monuclear reaction rates. 

Input files. Very often only the "input data file" (here- 
after IDF) is needed. It is read at the onset of the run and 
collects all the requirements needed for the calculations: 

— physical parameters: mass, chemical composition, mi- 
xing-length parameter, etc. 

— numerical parameters: maximum number of shells, 
kind of precision, etc. 

— criteria for halting the computations: age to be reached, 
value of the hydrogen abundance at centre, etc. 

— names and locations of the external data files con- 
taining the data of the tabulated EOS and opacity, 
names of the physical routines to be used, name of 
the model, of the set of units to be used, etc. 

^ Use of the PGPLOT package. 



The other input files have only specific functions, among 
them the most useful are: 

— "mixture" allows the use of an initial mixture not 
implemented in CESAM. 

— "modif jnix" allows to modify the abundances of some 
species in a mixture already implemented. 

— "reglages" allows to personalise the kind of preci- 
sion to be used (see below). 

— "planet" contains the characteristics of infall of plan- 
etoids. 

— "rap_iso" allows to modify the isotopic ratios. 

— "vent" defines the chemical composition of the wind 
when it differs from the atmosphere composition. 

— "zoom" allows to fit the resolution of the display for 
the plot on line. 

— "langue" allows to have the messages in English. 

The meaning of all items are explained in Paper 2. Ex- 
amples of files are given in the directory EXPLOIT. 

Output files. At the end of each time step, a "return 
binary file" (hereafter RBF) is created. It contains all the 
data needed to initialise or pursue a computation. There 
is the possibility, either to save all RBF, or to keep only 
the last RBF created. On request, "output data files" 
(hereafter ODF) are created. Three ODF are designed 
for adiabatic, non adiabatic and inversion asteroseismic 
investigations. These ODF also serve as input for some 
programs of the directory EXPLOIT. An ODF concerns 
the diffusion of the angular momentum (/3) . All items of 
output files are detailed in Paper 2. It is also possible to 
create personalised ODF. 

The kinds of precision. To optimise the calculations, sets 
of parameters, named "reglages", are fixed according to 
the kind of models to be calculated and to their subse- 
quent use. The most useful reglages are: 

— "realistic precision" for standard evolutions. 

— "super precision" used when a high level accuracy 
is needed. 

— "solar accuracy", close to super precision, but 
especially designed for seismological investigations; 
the number of shells of the last model is increased 
up to its maximal value. 

— "corot", close to super precision, but especially 
designed for investigations connected to CoRot. 

— "advanced" , for the computation of early type star 
models evolved up to the oxygen burning. 

— "normal precision" , for exploratory work. 

— "low mass", for the computation of late type star 
models. 

— "reglages", in that case, the parameters, designed 
by the user, are read on an input file named "reglages' 

The larger the expected accuracy, the larger the compu- 
tational expense. 
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Operating data and programs. The directory EXPLOIT 
contains: 

— examples of input files as IDF, "reglages", "planet", 
"modif _inix" , etc. 

— ASCII files of preliminary models for the initialisa- 
tions of PMS and ZAMS models. 

— miscellaneous programs to make plots of the chem- 
ical composition profile, or an extension of the grid 
on a given set of radius, or create the IDF for solar 
calibration, etc. 

The flow chart of CESAM is described in Sect. [2l As the 
numer ical features are detailed in the appendix of lMorell 
(|l997t PI (hereafter Paper I) they are only succinctly re- 
called in Sect. [31 except for the automatic allocation of 
mesh points described in Sect. 13.11 The restitution of the 
atmosphere is outlined in Sect. 31 The algorithms per- 
forming the temporal evolution are described in Sect. O 
The nuclear network is detailed in Sect. |6land the im- 
plementation of the rotation is described in Sect. [71 The 
various formalisms of convection implanted in CESAM are 
described in Sect. [9l Mass loss formalisms and infall of 
planetoids are described in Sect. [H EOS and opacities 
data available are listed in Sect. [TOl 



2 The flow chart. 

Initialisations. At the onset of the run, the IDF is read. 
Chemical composition is initialised according to the ini- 
tial mixture and to the isotopes used by the chosen nu- 
clear network. Then, using fit- formulas (see Sect. 16. ip . 
the thermonuclear reactions rates are tabulated on a rel- 
evant interval of temperatures0. Then, evolution begins: 

— either it starts from zero age on PMS or ZAMS: an 
initial model having the required specifications is de- 
duced from a model taken from a RBF or from a 
model in ASCII chosen in the directory EXPLOIT. 

— or it pursues a previous calculation, then the input 
is one RBF of the evolution going on. 

The evolution. The number of shells is updated as ex- 
plained in Sect. 13.11 Then, taking possible overshooting 
into account, the limits between radiative zones and con- 
vective mixed zones (hereafter LMR) are localised. The 
angular velocity (/3) and the chemical composition are 
then updated. In fine, the equations of quasi-static equi- 
librium for the interior and the atmosphere are solved. 
The process is repeated until convergence. 

When the nuclear engine is at work, the time step 
control is first based on the local accuracy achieved for 
the numerical integration of species of relevant interest 
(see Sect. I6.2l and Sect. I6.3|) and, second on a limitation 

^ Only available on electronic form. 

^ The errors introduced by these interpolations remain 
within the error bars of the data. 



of the relative changes of helium on the whole star. Oth- 
erwise, during the pre-main sequence, the limitation of 
the change of mechanical energy controls the time step. 
In case of divergence of any iterative algorithm the time 
step is halved. 

Stop criteria. According to flags, read in the IDF, the 
computations may be stopped: 

— when the expected age is reached. 

— when the central temperature reaches a given value. 

— as soon as the abundance of hydrogen at centre reaches 
a given value. 

— at the exhaustion of hydrogen at centre. 

— when the helium core reaches a given extent. 

— when the effective temperature crosses a given value. 

— at the ignition of the 3a cycle. 

— at the ignition of the carbon cycle. 

— at the ignition of the oxygen cycle. 

For most kinds of precision, the last time step is adjusted 
in order to fulfil the required stop condition. 



3 Numerical methods. 

Choice of variables. For the numerical integration of the 
stellar structure equations, the Lagrangian form is the 
most convenient as the discretisation on mass is readily 
expressed. However, it presents a singularity at centre 
and the core needs to be integrated apart. The Eulerian 
form of the equations does not suffer from such inconve- 
nience, but since the st ellar radius varies with t ime, there 
is a free boundary (see lStoer fc Bulirschl . ll979l par. 7.3). 
With the Lagrangian variables: ^ = M^/^ i?^, L^/^, the 
central singularity disappears (see [Paper J . par. Bl) and 
there is no need of a special treatment for the core (M 
is the mass, L the luminosity and R the radius). To be 
consistent, the chemical species are taken as functions 
of IVP^^. The pressure, P, and the temperature, T, are 
expressed in logarithms, (^ = InP, 77 = InT), on the 
ground that they change by more than six magnitudes 
from the centre to the atmosphere. 

Solving the differential equations. The unknowns are ap- 
proached by piecewise polynomials of order defined ac- 
cording to the required accuracy; the mostly used are of 
order 1, i.e. linear piecewise, and of order 2, i.e. parabolic 
piecewise. For the stellar modelling, such a flexible rep- 
resentation is well adapted to the presence of discontinu- 
ities resulting from the mixing of the convection zones. 
For the calculations, the piecewise polynomials are pro- 
j ected on a local linear basis of normalised B-splines 
(|De Booil . [19781 : [Schumakeil [l98lh . That allows to find 
back exactly the solution at any location. Moreover, B- 
splines basis are also used for solving: 

— the two points boundary initial value problems of the 
stellar structure and of the atmosphere by collocation 
(|De Booil [l978l ch. XV). 
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— the diffusion equations of c hemicals and anKula.r mo- 
mentum by finite elements (jOuarteroni fc Valiil . ll994h . 

The linear systems involved by the resolution of implicit 
equations are band-diagonal. 

However, due to the non-trivial and unfamiliar algebra 
of B-splines, the algorithms are much more elaborated 
than with the finite differences. Furthermore, efficient 
and stable algorithms have been constructed for inte- 
gration, differentiation, integration of differential equa- 
tions by collocation and for interpolation. In CESAM the 
routines, especially constructed to manage the calcula- 
tions with B -splines, are derived from the algorithms of 
ISchumake^ (1981, chap. 4). Details are given in .Paper I. 



3.1 Moving grid, mesh refinement and discontinuities 
tracking. 

An automatic mesh refinement is implemented. At time 
t, the locations of the n{t) mesh points are set by ful- 
filling the condition that, from a grid point to the next, 
the jump of a strictly monotonous "repartition function", 
Q(li,t), is equal to a "repartitio n constant" C(t) (see 
lEggletonL [l97l iPress et al.l . Il986l sect. 16.5). The loca- 
tions of grid points, /li, i = 1, . . . , n, known at the issue 
of the computations, satisfy: 



Q[ix,+i,t)~Q{ti,,t) = C{t), i^l, 



1. 



(1) 



The choice of (5(/i, t) is based on an a priori knowledge of 
the behaviour of the solution. For each t, one defines an 
"index" function q{p,t) mapping [0, /Xb] on [l,n]; the in- 
dex 1 {resp. n) corresponds to the centre and the index n 

2 

to the surface, i.e. = -^^cxt- Therefore, the integration 
is made on an equidistant grid. In terms of the derivative 
of Q with respect to g, Eq. ([T]) reads: 



The change of variables q{ii, t): 



9/i 
dq 



with e{^i, t) 



dQ 



is calculated from the analytic form of i). There are 
two more unknowns: ■(/'(t) and fi(q,t); they fulfil a sys- 
tem of differential equations of first order with boundary 
conditions: 



1, M = 





- Mb- 



The differential equations of internal structure and of 
atmosph ere, written with respect to q, are detailed in 
iPaper j . The equations are then solved on an equidistant 
grid, that allows numerical optimisations. 



Choice of Q. Q should be a strictly monotonous, two 
times differentiable function as simple as possible. By 
experiments, it has been found that the most convenient 
compromise is: 



Q^A^^ + Ar,r] + Af,^l, = A, = -1, A^ 
where ^ and rj have been defined above. 



15. 



Mesh refinement. The initial value of C{t) is fixed ac- 
cording to the expected level of accuracy. Along the evo- 
lution, its value is kept within ±2%, of its initial value 
by increasing (or decreasing) the total number of shells. 

Setting a grid point on a LMR. At the limit between a 
mixed zone and a radiative one, there may be a singu- 
larity of chemicals. Therefore each LMR has to coincide 
precisely with a grid point. To do that CESAM uses a 
weighted repartition function. On each side of a LMR, 
the weights are computed in such a way that they adjust 
locally the values of C{t) to the amounts just needed for 
iteratively "pushing" the closest grid point on the limit. 
In most cases, including mixing, the distances from the 
LMR locations to the closest grid points are lower than 
a few per cents of the characteristic local grid size. The 
more well defined the location of the LMR, the most 
efficient the algorithm. 

The grids. CESAM uses several grids for the B-splines 
representations of quasi-static variables, atmosphere and 
rotation variables. As it coincides with the adjustable grid 
of quasi-static variables, the Lagrangian grid used for 
the abundances of chemicals, is not fixed with respect 
to time. There is the possibility to use a "fixed grid", to 
avoid the numerical diffusion resulting from the variable 
Lagrangian grid. 



4 Atmosphere. 

The atmosphere connects the convective optically thick 
outer part of the envelope to the optically thin inter- 
stella r medium. In the interio r , the diffusion approxima- 
tion ( Kippenhahn fc Weigerti Il99l[ ) is used to simphfy 
the calculation of the radiative flux. In the outermost 
parts, this approximation is no longer valid, as soon as 
the Rosseland optical depth is lower than ~ 10: a special 
treatment is needed to restore the atmosphere. 
CESAM restores the atmosphere from a T{t,Tcs, g) law 
here, r is the Rosseland optical depth, Teff the effective 
temperature and g the gravity. The stellar radius, i?*, is 
defined as the bolometric one, i.e. the radius at the level 



where the local temperature is equal to Toff ([Morel et al.l . 

For genuine radiative T{t) laws, whatever is Teff, 
r* is located at a fixed optical depth, e.g. = 2/3 for 
the Eddington's law. For more precise laws, the value of 
is not fixed. Therefore, the location where Ti, is defined 
is a free boundary. 
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As oscillation modes are reflected in the outermost parts 
of stars, the pressure, the temperature and their gradi- 
ents should be continuous at the limit with the envelope. 
The continuity of the pressure gradient is trivially in- 
sured by the equation of quasi-static equilibrium verified 
on both sides of the limit. The continuity of the tem- 
perature gradient is more intricate to fulfil as the con- 
nection occurs in zones of convective instability. For the 
restoration of the atmosphere, the temperature gradient 
is derived from the T{t) law itself. Fixing the gravity the 
T(t) law is expressed as: = |t4j/(t). 



Using dr = —npdr and Viad = 

after some calculations: 
d\nT PkAT 



PkL 



~ dlnP 
3 PkL^, 



T dr GM 

2 



(2) 



R 



dr 



V, 



ad 



L 

IT 



R_ 



dr 



Here k is the Rosseland mean opacity, p the density, a the 
radiation constant, G the gravitational constant, the 
luminosity of the star, c the speed of light and Vrad the 
radiative gradient. In principle, the continuity is ensured 
when the same convection theory prevails on both sides. 
As in most cases the values differ, the continuity of the 
temperature gradient is insured by a weighted mean with 
respect to Inr. 

Eq. (12) is no longer valid for a genuine radiative T{t) 
law, as it ignores c onvection. In such a case, following 
the prescription of iHenvev et al.l (|l965l ) (M. Gabriel, J. 
Christensen-Dalsgaard, priv. comm.), the temperature 
gradient in the convective atmosphere is computed with 
a modified radiative gradient: 



For a genuine radiative T(t) law: 

lim f = 1, 

r>i dr 

therefore, at the limit between the atmosphere and the 
envelope, the radiative gradient is continuous and conse- 
quently the convective temperature gradient. 
The numerical integration of the differential equations 
fulfilled in the atmosphere is made by collocation, see 
iPaper j for details. The number of shells in the atmo- 
sphere is fixed from 50 to 100 according to the required 
level of accuracy. 



5 Evolution of the internal structure. 

Initial PMS model. The energy source in an initial PMS 
model is only of gravitational origin. At the onset of the 



Hayashi track, the star is fully convective, therefore isen- 
tropic and chemical ly homogen eous. The energy equa- 
tion is reduced as in llbeni (|1965D : 



dL 

Ml 



as 

'dt 



c{t)T 



(3) 



where c{t) is the "contraction constant" , and S the en- 
tropy. Along the interval of time At^ the energy radiated 
equals the change of gravitational energy, therefore, at 
first order: 



L{t) + L{t + At) 



At 



GAP 



At 



2GAP [R{t) 



R{t 
R{t^ 



-At) 
At)] 



GAP 



{L{t) + L{t + At))R{t)R{t + At) ' 



here R{t) is the stellar radius. An estimate of the initial 
time step is deduced from two models computed with 
Eq. ([2]), with close values of the contraction constant, 
c{t) a.TLilc{t+At). CESAM uses: c{t+At) = l.lxc(t). With 
c{t) = 0.02 cgs the temperature at centre is about lO^K; 
it is ten times larger with c{t) — 0.00008 cgs. Changes of 
the contraction constant allow to choose between differ- 
ent initial PMS models. Most of the PMS models can be 
initialised with a preliminary model in ASCII available 
in the directory EXPLOIT. CESAM assumes that a PMS 
model becomes a ZAMS model as soon as the release of 
gravitational energy balances the thermonuclear nuclear 
one. Such a model is chemically inhomogeneous. 

Initial homogeneous ZAMS model. A model of ZAMS 
with homogeneous chemical composition is not a physi- 
cal reality as the nuclear engine does not work at equilib- 
rium. However, it is a very convenient short way as, after 
a few time steps, the model is very close to the model 
at the end of the PMS. Several ZAMS initial models in 
ASCII are available in the directory EXPLOIT. 



6 Evolution of chemicals. 



6.1 The nuclear networks. 

At the onset of the computations, the abundances of 
chemicals are initialised according to the initial hydrogen 
and helium mass ratios and mixture. The initial abun- 
dance of each chemical species is split between its iso- 
topes, according to the isotopic ratios of nuclides. Sev- 
eral mixtures are impleme nted, the most useful are: the 
lAnders fc Greves sd fl98^ meteoritic mixture and the 



solar m ixtures of Greves se fc Noelsl ( 19931 ) and lGrevesse fc Sauvall 
(I1998D . - f necessary an IDF allows to modify the initial 
abundances of specific species or to use a mixture not 
yet implemented. 

Up to 16 nuclear networks are presently available. Hence, 
one can follow the evolution using only the chemical 
species and the thermonuclear reactions of interest. 
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The nuclear reaction rates are tabulated on relevant in- 
tervals of tem peratures. The rates are co mputed using 
the formulas of ICau ghlan fc Fowled (|l988f ) or of NACRE 
compilations (Angu lo et al.l. 1999 ) . For sol ar models, the 
improved rates of lAdelberger et al. ([l998^ are available. 
The weak screening of lSalpeted ( 196lh a nd the weak and 
intermediate screenings of iMitlen (|l997t) are available. 



6.2 Evolution without diffusion. 

The time scales involved in the temporal evolution of 
chemicals differ by a large number of magnitud es. From a 



here, V is the gradient operator and Vi the advection 
velocity. The symbol "•" means the vector inner product 

T 

and VM'-f 



dM 



The turbulent diffusion coefficients, are added to the i-th 
component of the vector, Di, of diffusion coefficients of 
the species i, see Eq. (fTH]) . 

For the integration, the abundances are approached by 
piecewise polynomials expressed on a B-spline basis with 
discontinuous derivat ives at each LMR . A fini te-elements 
method (see e.g. iQuarteroni fc: Vallil . 1 1994 ) is used to 
solve the diffusion equation. That allows an integration 
OTts which reduces to unity the order of the diffusion 



mathematical point of view it is a sfzff problem ( Hairer &: , , c^^x^tv^t. . „i i • r n • t •; mi i 
r-T-— I , , , , . ,1 1 ■ T~r equation. Ihe scheme is luily implicit, ihe nuclear term 

119911 ) and algorithms have been especiall y designed tor i , j c .^ ■ t i i r i 

T~T7T.^, ^ . . ..^ . Till I tt ■ c IS] evaluated as tor the implicit Eulcr s formula. 

it. Without microscopic diffusion. L-stable ( Hairer fc Wannea . . . j i ^ i i ^ i-o- • -i-i r 
■777^ , „N . ., 1 ^ 7, The mixing is made by turbulent diffusion with coet- 



I1991L par. 4.3) implicit Rungc-Kutta schemes, are avail- 
able for the chemical evolution. 

In the radiative zones, the equations to be solved are 
formally written: 

OX' 

-^ = <F,iP,T,X,t), l<i<nx. (4) 

Here, Xi is the abundance per mole of the chemical species 
i, Wi the rate of change of a;^, X the vector of chemical 
abundances, t the time and rtx the number of chemicals; 
Xi is given by: 
_ 

Xi is the abundance mass fraction and Vi the atomic 
mass. In mixed zones (hereafter MZ), convective eddies 
homogenise the chemical composition. There, the mix- 
ing and updating of chemicals are done simultaneously, 
therefore the changes of mean abundances read: 

^,{P,T,X,t)dM / [ dM. (5) 
™ Jmz I Jmz 

As a grid point is defined on each LMR (see Sect. 13. ip . 
the discontinuities of the abundances are explicitly cal- 
culated. 

Control of the accuracy. A good estimate of the nu- 
merical accuracy o f an integration is obtai ned with the 
Fehlberg method (|Stoer fc Bulirschl . Il979l par. 7.5.2). 
It needs to triple the calculations. As it is prohibitive, 
falling anything better, the time step is simply adjusted 
in such a way that, over a time step, the relative changes 
of the abundances remain within fixed limits. The largest 
the expected accuracy, the narrower the limits.. 



6.3 Evolution with diffusion. 

With microscopic diffusion, the equations of the evolu- 
tion of chemicals have the form: 



dxj 
dt 
F, 



m 

dM 



dxi 
'dt 



(6) 



nucl. 



ficient Dmz » 1- At each LMR, the abundances Xi 
and fluxes Fi are continuous functions with discontinu- 
ous first derivatives, owing to the jumps of the diffusion 
coefficients. Therefore Eq. © holds everywhere. Two 
formalisms are available for the calculation of the diffu- 
sion vector Di'. 

— The co efficients are calculated according to lMichaud fc Proffittj 
(|1993[) . The metais are "test elements" , their diffusion 

only results from collisions against protons. Based on 
the presence of protons, this formalism is only valid 
for the main sequence. 

— The diffusion c oefficients are computed according to 
iBurgersI ( 1969t) . this formalism is outlined beneath. 

Boundary conditions. At the outermost limit, Al — Afext, 
it is assumed that there is neither input nor output of 
matter, then Fi(Mext) = for any particle i. At centre 
Al = 0, because of the spherical symmetry, Fi{Q) = 0. 

Control of the gravitational settling. For stellar models 
with mass larger than w 1.4 Mq, the use of microscopic 
diffusion alone produces an important depletion of he- 
lium and metals at the surface and a concomitant en- 
hancement of hydrogen. Different ways to overcom e this 
problem have been used. lEggenberger e t al. ' ('2005') intro- 
duce some turbulence due to rotation. Di Mauro J200^ 
suppr ess diffusion in the outer layers. IChabover et all 
(|1999l ) include a win d mass loss whi c h red uces the dif- 
fusion in outer layers, iTurcotte et al. 1 (|1998D introduce a 
turbulent mixing. CESAM allows to control gravitational 
settling with a radiative turbulence of coefficient d^, pro- 
portional to the radiative kinetic viscosity; it results from 
the energy exchanges between thermal collisions leading 
to excitation and ionisation of atoms an d ions (Thoma_^, 
I1930I: iMihalas fc Weibel-Mihalad . [T983 . p. 461-472). 



Re^- 



4 aT^ 



(7) 



15 CKp^ ' 

c is the speed of light in vacuum. The phenomenolog- 
ical parameter Rei, has been found close to unity by 
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iMorel fc Thevenin I (|2002h . The physical meaning of this, 
as efficient as s imple source of turbulent m ixing, has been 
questioned by lAlecian fc MichaudI HqoI). 



Burgers 's flow equations. With respect to the abscissa x, 
the density number ni and Wj, the diffusion vel o city o f 
the particle i, are related bv llben fc MacDonaldl ( 19851 ): 



'dt 



dx 



'dt 



nucl. 



With the formalism of iBurger^ ( 1969[ ). the diffusion ve- 
locity is expressed as: 



Kij are the so-called resistance coefficients, they repre- 
sent the effects of collisio ns between particles i and j 
(iMichaud fc ProffittL[l993h : 



16 



(11) 



(12) 



with rriij = {mimj)/{mi + rrij) as the reduced mass of 
particles i and j. The heat flux terms involve additional 
collision integrals: 



2 n 

= 1 — _ ^-^ 



(12) 



5 o(") 



2.5 



2 5a 



(12) 



(13) 



a 



(22) 



a 



(11) 



a 



(11)' 



where the diffusion coefficient of the particle i with re- 
spect to the particle j, bij, and the advection velocity, 
Vi, come from the solution of a linear s ystem. The dif- 
fusion velocities , for ions and electr ons, ( Burgerd . Il969l : 
I Cox et aT ], 119891 : iThoul et all . \l994 , satis^H 



dR 



P^dP 
p dR 



— — riiZjCh = 



5 , dT 
2"^^di? 



i 



j 

m^r. 



Kij{wj - Wi) + 
rrnrj 



rrii 4- wij 
mj[wj ~ Wi) 



(8) 



— K 



r, V K, 



3mf -f m]z[^ -I- \m,mjz'l^ 



5 v 



{rrii 



(9) 



e is the electron charge, E the electric field, Pi the par- 
tial pressure, rui — Vivriu is the mass of particle z, k the 
Boltzmann constant, and the magnitude of the resid- 
ual heat flow vector; pi the partial density is given by: 



Pi = pXi = pXiUi = niV^rUu, 



(10) 



with TVq as the Avogadro number, the inverse of the 
atomic mass unit m„: 



X,, 



Ui = pNo — = pNoXi. 



(11) 



The force equation ^ represents the pressure and the 
concentration dependence of the diffusion velocity, while 
Eq. ([9]) prevails for its thermal dependence. The charge 
of any isotope is taken as the averagecjl charge, zi , over 
all its ionisation states. An unique mean charge is used 
for all the isotopes of a given chemical. The quantities 

* Without conduction current and magnetic field. 
^ Weighted by the ionisation rates. 



iPaquette et al.l (|1986( ) showed that a|j^'^ can be written: 
aif)^^;f)e,„ (fcO-(ll), (12), (13), (22), 



e 

T 



For attractive and repulsive screened Coulomb poten- 
tials, the quant ities ln(F'^'^'^)y have been tabulated by 
iPaquette et al.l ()1986) . The equations of dynamical con- 
servation of the mass and charges respectively are: 



nx 



i=l 



XiViWi 



nx+1 
0, ^ ZiXiWi 

4=1 



0, 



(13) 



where nx + 1 is the index of el ectrons. Following Paauette et al.l 
(I1986D and other works, e.g. | lben fc MacDonaldl (11985^ 
ICox et al.l (|l989t) . iThoul et ah (T994.) . the diffusion ve- 
locities Wi come from the solution of the system of 2nx -\- 
3 linear equations, formed by the 2nx Eq. ([8]) and Eq. ^ 
for the ions, the Eq. ([9]) for the electrons and the two 
Eq. (|13p . The 2nx + 'i unknowns are Wi, i = 1, . . . , nx + 
1, Ti, i — 1, . . . , nx + 1 and E. 

For want of something better, with the ideal gas law, the 
pressure and the partial pressures respectively are: 

p-RT 



P 



Pi = nikT = pxiNokT = Ppxi, 



(14) 



here TZ is the perfect gas constant, and p the mean molec- 
ular weight. In spherical symmetry, the pressure and the 
temperature gradients are given by: 



dP 
dR 



dT T dlnT dV 



-P'^ dR 



(15) 



V dhiV dR 

Working with Eq. (fTO|) to Eq. (fT5|) , with respect to p and 



Xi, the Burgers's equations (p 
,dxi 



may be rewritten: 



dXn 



dR' 



OR 



,0,...,0)^ 



The solutions are: 
Lu^V + BD^, V 



B = A'-^G. 
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For abridgment, neither the derivation of above equa- 
tions, nor the comphcated forms of vector 7 and matrix 
A and G, are reproduced, all details are given in Paper 3. 
The diffusion velocities of the ions are expressed as: 

nx 

Wi = '^hj-^ + Vi, i ^ 1, . . . ,nx, 
i=i 

and, owing to Eq. (fTTjl . the diffusion vector reads: 

A = {xibii, . .., Xibinj, )^ , (16) 

for i, j = 1, . . . , nx, bij and w^, are respectively the coef- 
ficients of matrix B and vector V. 



Table 1 Comparison between the mean charges from the 
approximate solution and the precise equation of state of 
iGabriel (:1997i ) (label "G") beneath the convection zone and 
at centre, for a calibrated solar model, . 



T=2.2 MK centre 
elements Zg Z Zq Z 



G 


5.89 


5.92 


6.00 


6.00 


N 


6.75 


6.82 


7.00 


7.00 





7.47 


7.56 


8.00 


8.00 


Fe 


17.0 


16.7 


24.2 


24.2 



Calculation of mean charges. The Burgers formalism in- 
volves the charges of the isotopes. To simplify, CESAM 
considers a unique mean charge for all the isotopes of 
each chemical. For the calculat ion of the ion i sation rates, 
the Saha-Boltzmann equation ( Cox fc Giulil . ll968l . eq.l5- 



30) has been adapted in the following way. Let rij^i be 
the number density of atoms in ionisation state j of the 
chemical species i. The ratio of the total number of atoms 
in successive stages of ionisation can be written: 



the Qj^i are the statistical weights of fundamentals. For 
the smoothing function /(x), CESAM uses the piecewise 
cubic polynomial with zero derivatives at a; = and 
X — p {p — A): 



nj-i.t 



■ exp 



(17) 



here Xj,i is the ionisation potential, Uj^i the partition 
function and rj the electron degeneracy, related to the 
number density of free electrons and te mperature, thr ough 
the ha lf integer Fermi-Dirac function ( Clavtonll 19681 eq. 
2-57). lEggleton et al.l ( 19731 ) have introduced a conve- 
nient approximate treatment of pressure ionisation: a 
numerical correction is postulated to ensure that the 
plasma remains completely ionised at sufficiently high 
density. In the inner solar radiative zone, this approxi- 
mate treatment leads to too large mean charges for the 
ions and, at the centre, iron is f ully ioni s ed, w hile it is 
only 85% according to Table 1 in lGabriell (|l997l ). A sim- 
ilar behaviour is observed with the modifie d parameters 
recommended bv lProffit fc MichaudI (ll99lL eq. 4). 
In CESAM, the partition functions are limited to the sta- 
tistical weights of fundamental levels and, as soon as the 
mean distance between the ions becomes of the order 
of the size of the ion cloud, i?D, the ratio of statistical 
weights is sm oothly reduced to zero. The Debye-Huckel 
radius writes (|Clavtonl . fl968l eq. 2-235) *: 



i?n = 



kT 



The Saha-Boltzmann Eq. (fTTl) is then written: 
x=^~l, x'j.i = max(0, x»j - xfj), 



f{x)e^^(r,+ ^) 



(18) 



c ^ Je 







fix) 



Then, 

















+ 3 











1 



if a; e [0,p] 
if X > p, 



as soon as x^j < X?j, the quantity gj^i,i/ gj^J{x) 
goes to zero and the level j — 1 becomes fully ionised. 
The set o f Saha-Boltzma nn equations for all ions is writ- 
ten as in iMihalasI ([1978, eq. 5-17). They are solved by 
iterations using a second order Newton-Raphson scheme. 
For the mean charges, quantities of main interest in the 
present investigation. Table [T] reveals agreements bet ter 
than ±3% with the data of Table 1 of lGabriell (|l997l ). 



7 Rotation. 

Rotation is considered in CESAM under the assumption of 
spherical symmetry. With non zero angular velocity, the 
mean centrifugal acceleration affects the local gravity. In 
the initial model, rigid rotation is assumed. The initial 
angular velocity can be read from the IDF in different 
units: 

— radian/s. 

— km/s, it corresponds to the rotational velocity of the 
outer part of the star. As the outer radius depends 
on the outer gravity, the initial model has to be iter- 
atively adjusted to have the required outer velocity. 

— days, that corresponds to the rotation period. 



7.1 Rotation without diffusion of angular momentum. 
Several options are available: 

— no rotation: the initial angular velocity must also be 
zero. 

— solid-body rotation where the angular velocity is kept 
to its initial value. 
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— rotation is solid and angular momentum is globally 
conserved. The angular velocity changes with respect 
to time, according to structural changes. 

— angular velocity changes according to the conserva- 
tion of the local angular momentum. The rotation is 
not solid, but in the mixed zones. 



border of a convective core, the discontinuities of chem- 
icals are assumed to be eroded only by diffusion. With- 
out diffusion (/3), the abundances profiles are spatially 
linearly interpolated between their values on the convec- 
tive core and at the former location of the LMR. That 
avoids a noisy behaviour of chemical gradients and, con- 
sequently, of the profile of the Brunt- Vaissala frequency. 



7.2 Rotation with diffusion of angular momentum. 



I I , I ID Equation of state and opacity. 

The tw o formalisms of lTalon et all (|l997l ) and lMathis k. ZahnI 

of the diffusion of angular momentum are im- 
plemented in CESAM. The diffusion coefficients of the 
angular momentum are c omputed acco r ding either to 
iPalacios et all (|2003( ) or to lMathis et"all (|2004l ). 



8 Mass loss and infall of planetoids. 

Several formalisms of mass loss are implemented. The 
mass loss rate is either negative or positive, i.e. increase 
of mass. With diffvision, the chemical composition of the 
input [resp. output) can differ from those of the outer 
convective zone, from where the output {resp. input) is 
assumed to come from {resp. vanish). The mass loss rate, 
in unit of MQ.yr~^ is read in the IDF. The following 
options are implemented: 

— standard mass loss. 

— solar mass loss, the mass loss is halted as soon as the 
mass of the model reaches the solar value. 

— mass changes due to nuclear energy generation. 

— infall of planetoids (/3). The characteristics of the in- 
fall, duration, amount of terrestrial masses, chemical 
composition of the planetoids, are read in a IDF. 

Angular momentum losses are implemented but in vali- 
dation {(3). 



Four an alytical EOS are impl emented. The most useful 

are e ff (|Eggleton et al.l . [T973[ ) and CEFF ([Christensen-Dalsgaard fc Dapt 
'1992^). Numeri cal EOS are available. The mhd tables 
(Mihalas et al.l il9881 and the opal 1993 and opal 2001 
tables are used with the OPAL interpolation scheme for 
tables with Z = 0.01 and 0.02. 

OPAL 1996 opacity tables are implemented, with Kurucz 
low-temperature values. The ratios between the abun- 
dances of heavy-elements are fixed at their initial values, 
regardless of changes due to nuclear reactions and dif- 
fusion. All OPAL tables of type 2, available today, are 
implemented altogether with z14xcotrin21, the pack- 
age of A.I. Boothroyd. 

The neutrinos emitted in nuclear reactions are taken 
into account, it is assumed that they can freely escape 
the star. Other physical processes related to neutrino 
such as Urea process, plasma neutrino, pair neutrino 
(jKippenhahn &: Weigertl . [l99ll par. 18.6) have only been 



implemented in private codes derived from CESAM. 
For temperature values above T > 7 lO^K ~ 7Kev, the 
plasma is fully ionised, so the Rosseland mean opacity 
is reduced to the Compton scattering by free electrons 
(Cox fc aiulL,1968. par. 16.6). 



11 Calibration of the Solar model. 



9 Convection. 

Two formalisms for the computation of the temperature 
grad ient in the convection zones are available: the stan- 
dard [Bdh^^Vitensi il95&) mixing-length (MLT) formal- 
ism is considered with the optical thickness of the con- 
vective bubble, and the lCanuto fc Mazzitellil (|l99ll ) for- 
malism. The mixing-length parameter is read in the IDF. 
Overshooting beneath and/or above the convection zones 
can be accounted for. The overshooting parameters, scaled 
by the local pressure scale height, are read in the IDF. 
In overshoot regions, the temperature gradient is set, ei- 
ther to the adiabatic or to the radiative gradient. The 
convective zones and their extents by overshooting are 
homogenised, see Sect. 16.21 

Up to now specialised treatment of the semi-convection is 
not implemented in CESAM. With microscopic diffusion, 
in areas swept across by the backward movement of the 



We have performed calibrations of the solar model with 
CESAm2k with various sets of input physics and ini- 
tial parameters. This consists in adjusting the initial pa- 
rameters of the model (helium abundance and mixing- 
length parameter) in order to satisfy the observational 
constraints on the solar global parameters at solar age. 

We have adopted the values of the astronomical and 
physical constants specified for the calculation of the stel- 
lar models com pared in the different ESTA tasks (see 
iLebreton et al.L ,2007). For the solar global parameters, 
we therefore took Rq = 6.9599 10^° cm (radius), Lq — 
3.846 10^3 erg.s-i (luminosity) and Mq = 1.98919 lO^^ g 
(mass) . The value Rq refers to the radius of the model 
layer where T = Tcs = 57 77 K. For solar age, we a dopted 
the value = 4.57 Gyr (|Bahcall fc et al.l . ll995D . 

Table [2] presents the sets of input physics used in the 
different calibrations of solar models. Al l models have 
been calculated with the OPAL 2001 EOS (.Rogers fc Navfonovl . 
and the 1995 OPAL interior opacities (jlglesias fc Rogersl 
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Table 2 The sets of parameters used for the calibration of solar models with CESAM2k. We used (i) two solar heavy 
elements mixtures, the GN93 (jGrcvcssc & Noels, 1993) or the AGS05 (Asplund ct al., 2005) mixture; (ii) 1;wo formalisms for 
convection treatment, the BV (Bohm-Vitcnsc, 1958) and CGM (jCanuto ct al., 1996) formalism, (iii) two T(r) laws for the 
atmosphere, the Eddington grey law and T(r) laws derived from ATLA S9 model atmospheres (Kurucz, 1992) calculated with 
the same convection formalism as the interior models ( seelSamadi et al.. 2006.') a nd (iv) tw o formalisms fo r the calculation 
of the microscopic diffusion of the elements, the MP93 (jMichaud fc ProfHttl . Il993l ) and B69 (|Burgersl . [19691 ) formalism. 



Model 


Mixture 


Convection 


T(r) law 


Diffusion 


A 


GN93 


BV 


ATLAS9 


MP93 


B 


GN93 


BV 


EDDINGTON 


MP93 


C 


GN93 


BV 


ATLAS9 


B69 


D 


GN93 


CGM 


ATLAS9 


MP93 


E 


AGS05 


BV 


EDDINGTON 


MP93 



Table 3 Properties of the five CESAM2k calibrated solar models. Yq, Zq and {Z/X)o are the initial values of the helium 
and heavy elements mass fractions and (Z/X) ratio. Yc, Za are the present values in the convective envelope while i?c is the 
radius at the basis of the convective envelope in units of the solar radius, pc, Tc and are the central values of the density 
in g.cm~^, temperature in 10® K and hydrogen abundance. Qconv.int is the value of the convection parameter in the interior 
(q = I /Hp where I is the mixing length and Hp is the pressure scale-height) 



Model 


Yo 


Zo 






z. 


Rc 


Pc 


Tc 




Q^conv.int 


A 


0.2735 


0.0196 


0.0278 


0.2447 


0.0181 


0.7143 


153.1 


15.72 


0.3387 


2.42 


B 


0.2735 


0.0196 


0.0278 


0.2447 


0.0181 


0.7145 


153.1 


15.72 


0.3387 


1.76 


C 


0.2741 


0.0198 


0.0280 


0.2460 


0.0180 


0.7152 


153.0 


15.72 


0.3388 


2.40 


D 


0.2735 


0.0196 


0.0278 


0.2447 


0.0181 


0.7143 


153.1 


15.72 


0.3387 


0.77 


E 


0.2637 


0.0141 


0.0196 


0.2339 


0.0129 


0.7297 


151.0 


15.52 


0.3578 


1.72 



0.015 



0.01 



0.005 



-0.005 




0.2 



0.4 



0.6 



0.8 



Figure 1 Rela tive differences be tween the seismic sound 
speed derived bv lBasu et al.l l|2000l ) and the solar models pre- 
sented in Table [21 



Il996l) . T wo sets of low temperat u re op acities have been 
use d: the Alexander fc FergusonI ( 1994[ ) tables giv en for 



erations) but we considered eith er the Michaud fc ProffifrO 
(|1993| . hereafter MP93) or the |Burged ilMSl, hereafter 
B69) formalism. Co nvection is treated ei ther according 
to the cla ssical MLT (iBohm- Vitensd . [Tosl hereafter BV) 
or to the ICanuto et al. ( 19961 . hereafter CGM) formal- 
ism. In the CGM fo rmalism, like in the CM formalism 
(jCanuto fc Mazzitelli . J,991), the contribution of eddies 
with different sizes is taken into account in the calcu- 
lation of the conve ctive flux and velocity. In addition 
ICanuto et al. I (ll996i) take into account the feedback of 
the turbulence on the energy input from the source which 
generates turbulent convection. For the atmosphere cal- 
culation we considered either the classical Eddington 
grey r(r)-law or a law d erived from Ku rucz's ATLAS 9 
1-D model atmospheres ( Kurucd . Il992!) . We have taken 
the sa me r(T)-laws as used in the work bv lSamadi et al.l 
(120061 ). These laws are based on model atmospheres cal- 
culated with either the BV or the CGM convection for- 
mulation. In both cases, the atmosphere calculation was 
performed adopting a value of the mixing-length param- 
eter aconv.atm = 0.5 which allows to fit a t best the ob- 

serve d profiles of the solar Balmer lines (see lvan't Veer-Menneret and Meg 
Il996l ). Therefore aconv.atm is different from the value of 
aconv.int in the interior, this latter being adjusted to cal- 
ibrate the solar model. Finally, we adopted the GN93 so- 
lar mixture of heavy elements (.Grevesse fc Noels, 1993) 



thelGrevesse fc Noelsl (|l993l) solar mixture and the lFergusoftnefe 4Msplund et all . [20051 ) which is derived from a time- 
(|2005^ ta bles g i ven fo r the new solar mixture derived by dependent, 3-D hydrodynamical model of the solar atmo- 
lAsplund et a L ('2005') (see below). All models take into sphere. The abundances of C, N, O of the AGS05 mix- 
account the diffusion of chemicals due to pressure, tern- ture are smaller than in the GN93 one which leads to an 
perature and concentration gradients (no radiative accel- 
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important decrease of the solar {Z / X) ratio: {Z/X)q ~ 
0.0245 for the GN93 mixture and {Z/X)q = 0.0171 for 
the AGS05 mixture. 

Table [3] presents the results of the solar model cal- 
ibrations. The relative differences in radius, luminosity 
and present [Z/X) surface value of the five models A, 
B, C, D, E with the observed values are lower than 
10~^. The relativ e differences between the seismic sound 
speed derived by iBasu et all 12000.) and the models are 
plotted in Fig. [TJ The helioseismically measured values 
of the present radius at the base of the convective en- 
velope i?e,o and of the present solar envelope helium 



abundanc e Y^ o provide 



strong constraints for the so- 



lar model. iBasu fc Antia (Il997) hehoseismically derived 
i?e,o = 0.713±0.001i?g. lBoothrovd fc SackmannI jlool) 
derived a mean value 1^ o = 0.245 ± 0.005 from differ- 
ent helioseismic determinations. In all our models but 
one (model E), we find values of Yg and Re in reason- 
able agreement with the seismic values. Model E is based 
on the AGS05 solar mixture which makes the agree- 
ment between the solar model and helioseismic observa - 
tions much worse (see for instance [Basu fc Antial . [200^ . 
More details on the solar models calculated with CESAM 



and th eir seismic prop erties can be found in Morel et al.l 
(|l999D : lProvost et all (|2000l ); IZaatri et al.l (|2007^ 
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